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In the preceding paper ( cond-mat/0405252 ), we have conjectured that the main transport prop- 
erties of a dilute gas of inelastic hard spheres (IHS) can be satisfactorily captured by an equivalent 
gas of elastic hard spheres (EHS), provided that the latter are under the action of an effective drag 
force and their collision rate is reduced by a factor (1 + a)/2 (where a is the constant coefficient of 
normal restitution). In this paper we test the above expectation in a paradigmatic nonequilibrium 
state, namely the simple or uniform shear flow, by performing Monte Carlo computer simulations of 
the Boltzmann equation for both classes of dissipative gases with a dissipation range 0.5 < a < 0.95 
and two values of the imposed shear rate a. It is observed that the evolution toward the steady state 
proceeds in two stages: a short kinetic stage (strongly dependent on the initial preparation of the 
system) followed by a slower hydrodynamic regime that becomes increasingly less dependent on the 
initial state. Once conveniently scaled, the intrinsic quantities in the hydrodynamic regime depend 
on time, at a given value of a, only through the reduced shear rate a*(t) oc a/y/T(t), until a steady 
state, independent of the imposed shear rate and of the initial preparation, is reached. The dis- 
tortion of the steady-state velocity distribution from the local equilibrium state is measured by the 
shear stress, the normal stress differences, the cooling rate, the fourth and sixth cumulants, and the 
shape of the distribution itself. In particular, the simulation results seem to be consistent with an 
exponential overpopulation of the high-velocity tail. These properties are common to both the IHS 
and EHS systems. In addition, the EHS results are in general hardly distinguishable from the IHS 
ones if a > 0.7, so that the distinct signature of the IHS gas (higher anisotropy and overpopulation) 
only manifests itself at relatively high dissipations. 

PACS numbers: 45.70.Mg, 05.20.Dd, 05.60.-k, 51.10.+y 



I. INTRODUCTION 



A granular gas in the rapid flow regime is usually mod- 
eled as a system of smooth inelastic hard spheres with a 
constant coefficient of normal restitution a. The key in- 
gredient of this minimal model is that energy is not con- 
served in collisions: in every binary collision, an amount 
of energy proportional to 1 — a 2 is transferred to the 
internal degrees of freedom and thus it is lost as trans- 
lational kinetic energy. Therefore, the gas "cools" down 
and the granular temperature monotonically decreases in 
time unless energy is externally injected into the system 
to compensate for the collisional loss. If this energy injec- 
tion takes place through the boundaries (e.g., vibrating 
walls, thermal walls, sliding walls, . . . ) a nonequilibrium 
steady state can be reached characterized by strong spa- 
tial gradients in basic average quantities, such as density, 
mean velocity, or temperature. Kinetic theory tools can 
be straightforwardly extended to granular gases and, in 
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particular, the Boltzmann and Enskog kinetic equations 
have been formulated for inelastic collisions Q, 0, ■ 

In the preceding paper 0, we have suggested that 
the nonequilibrium transport properties of a genuine gas 
of inelastic hard spheres (IHS) can be accounted for, 
to some extent, by an "equivalent" gas of elastic hard 
spheres (EHS). This requires the inclusion of two basic 
points. First, the EHS are assumed to feel the action of a 
drag force with a friction coefficient that mimics the col- 
lisional cooling rate of the true IHS gas. This guarantees 
that the local energy balance is approximately the same 
in both systems. Second, the collision rate of the EHS gas 
must be decreased by a factor f3(a) with respect to that of 
the IHS gas at the same (local) density and temperature. 
This can be interpreted under the assumption that, while 
the EHS have the same mass m as the IHS, they have a 
diameter a 1 = /3 1//2 c (where henceforth we are restricting 
ourselves to three-dimensional systems) smaller than the 
diameter a of the IHS. Comparison between some basic 
collision integrals for IHS and EHS suggests the simple 
choice (3 = |(1 + a) Q to optimize the agreement be- 
tween both descriptions. 

The aim of this paper is two-fold. First, we want to 
assess to what extent the EHS gas succeeds in mimicking 
the transport properties and other quantities of the IHS 
gas by choosing the most widely studied nonequilibrium 
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state, namely the simple or uniform shear flow (USF). 
Since the USF state is intrinsically non- Newtonian 0, , 
it provides an interesting test case to check whether or 
not the IHS and EHS gases behave similarly in extreme 
situations far from equilibrium. As we will see, it turns 
out that most of the noncquilibrium properties of both 
types of system are practically indistinguishable for de- 
grees of dissipation as large asa~ 0.7. For higher dissi- 
pations, the agreement is still at least semi-quantitative. 
The second goal of the paper is to carry out a rather ex- 
tensive study of the most relevant properties of the USF 
in dissipative gases. This study includes the unsteady or 
transient regime (usually neglected in previous studies), 
which can be decomposed into a kinetic stage followed 
by a hydrodynamic stage. In the steady state we pay 
special attention not only to the rheological properties, 
but also to the velocity cumulants and to the distribution 
function itself. 

The paper is organized as follows. The Boltzmann 
equation for both types of system, IHS and EHS, is briefly 
recalled in Sec. |n] where also part of the notation is in- 
troduced. This is followed in Sec. HHI bv a description of 
the USF and of the main quantities of interest. Section 
II V I is devoted to some details of the numerical simulation 
method employed to solve the Boltzmann equation for 
each class of systems. The most important part of the 
paper is contained in Sec. [V] where the results for the 
transient and steady-state problems are presented and 
discussed. When possible, the simulation data are also 
compared with predictions from a simple kinetic model 
• The paper is closed by some concluding remarks 
in Sec. |VT| 



II. BOLTZMANN EQUATION 

In a kinetic theory description the relevant quantity is 
the one-particle velocity distribution function /(r, v;i). 
Its first five moments define the number density 



n(r,t) = / dv/(r,v;t), (2.1) 

the nonequilibrium flow velocity 

u(r,t) = /dvv/(r,v;t), (2.2) 
n(r,t) J 

and the granular temperature 



T(r,t) 



m 



3n(r, t) 



dvU 2 (r,<)/(r,v;t), (2.3) 



where m is the mass of a particle and V(r, t) = v — u(r, t) 
is the peculiar velocity. The hydrostatic pressure p = nT 
is one third the trace of the pressure tensor P defined as 

P(r,t)=m J dvV(r,i)V(r,t)/(r,v;i). (2.4) 



A. Inelastic hard spheres 

The evolution of / in the low-density limit is governed 
by the Boltzmann equation. In the case of a gas of in- 
elastic hard spheres (IHS), it reads 0,0,01 

(9 t +vV)/ = jW[/,/], (2.5) 
where [/, /] is the Boltzmann collision operator 

x[a- 2 /(v")/K)-/(v)/( Vl )], (2.6) 

the explicit dependence of / on r and t having been omit- 
ted. In Eq. (|2.ti|l . a is the diameter of a sphere, is 
the Heaviside step function, <x is a unit vector directed 
along the centers of the two colliding spheres at contact, 
g = v — vi is the relative velocity, and a is the coefficient 
of normal restitution. The pre-collisional or restituting 
velocities v" and v" are given by 

v"=v -_ (g-o-)o-, vi'=vi + - SI -(g-ff)o- > (2.7) 



2a 

while the direct collision rule is 
l + a 



2a 



v = v— ■ 



l + a 

(g-CT)a-, v x =ViH — (g-ff)a-. (2. 



The inelasticity of collisions contributes to a decrease 
of the granular temperature T(t), i.e., 



m 

3/7 



J dvV 2 jW[f,f]=-(T, 



(2.9) 



where £ is the cooling rate. By standard manipulations 
of the collision operator, the cooling rate can be written 
as 00 



(2.10) 



where 



(V&) = ^J dv iJ dv 2 l v i - v 2 | 3 /(vi)/(v 2 ) (2.11) 

is the (local) average value of the cube of the relative 
speed and 



A 



v / 2T/m 
is a (local) characteristic time, 

(V2 nnc 



X 



(2.12) 



(2.13) 



being the (local) mean free path. Note that the cooling 
rate £ is a nonlinear functional of the distribution func- 
tion / through the average (Vfy) an d cannot be explicitly 
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evaluated without the knowledge of /. Nevertheless, a 
simple estimate is obtained from Eq. I|2.11[l by replacing 
the actual distribution function / by the local equilibrium 
distribution 

/o(v) = n(m/27rT) 3/2 exp(-mF 2 /2T). (2.14) 
In that case, 

3/2 



32 T 



(2.15) 



Insertion of this approximation into Eq. (|2.1UI) yields the 
local equilibrium cooling rate 0, @ 



Co 



2t- 



3^ 



(1-a 2 ). 



(2.16) 



This local equilibrium estimate is still a functional of /, 
but only through the local density and temperature, i.e., 
Co oc nT 1 / 2 . In addition, its dependence on inelasticity 
is simply ( ocl-a 2 . 

The characteristic time r defined by Eqs. (|2.12(1 and 
(12.13(1 is of the order of the (local equilibrium) mean free 
time T m ft, namely 



A 



T m ft 



<y)t 



-T, 



(2.17) 



where in the last step we have taken into account that 
the mean (peculiar) speed is (V) — > (V)o — y/ST/nm in 
the local equilibrium approximation. It is also convenient 
to introduce a characteristic time t v associated with the 
momentum transfer or viscosity. Its expression is 



r v = — = 1.016— r, 



(2.18) 



where r/ = 1.016 x 5 yj mT / V / 16a 2 is the Navier-Stokes 
shear viscosity in the elastic limit (a — > 1) |9j. Note that 
t v /t f» r/r mft w 1.13. 



B. (Frictional) elastic hard spheres 

Now we consider a dilute gas of elastic hard spheres 
(EHS) of the same mass m as the IHS but with a diam- 
eter a' = [/3(a)] "^cr smaller than that of the IHS. As a 
consequence, the characteristic time r' and the mean free 
path A' of the EHS gas are [l(j 



A' 



nira 



(2.19) 



Furthermore, we assume that the EHS are under the in- 
fluence of a drag force Fd rag = — TO7(a)V with a friction 
constant 7(a) = ^Co(a), where Co(o0 is given by Eq. 
(|2.16() . Therefore, the Boltzmann equation for the (fric- 
tional) EHS gas is 

U + v . V - ||- • v) / = (3 J« [/, /], (2.20) 



where the elastic collision operator [/, /] is given by 
Eq. I|2.6I) by setting a — 1 both explicitly and in the col- 
lision rule l|2.7[) . but keeping the factor a 2 . Henceforth, 
when referring to the EHS system, we will always under- 
stand that the particles are frictional, in the sense that 
the external force Fdrag is acting on them ^l| • 

Friction produces in the EHS gas a cooling effect char- 
acterized by the rate (,o{a). This is intended to mimic 
the cooling effect in the IHS gas due to the collisional 
inelasticity, which is characterized by the rate C(a) given 
by Eq. I|2.10|l . Both cooling rates are quantitatively close 
each other as long as the density and temperature are 
similar in both systems and l|2.15[) is a good approxi- 
mation. In principle, we could consider a friction con- 
stant 7 = |Co(^i2)/(^i2)o) so that the cooling rate of 
the EHS would be the same functional of / as in the 
IHS case. However, this would complicate excessively 
the EHS model without a correlated gain in accuracy, as 
we will see in Sec. IV CI 

The parameter (3{a) = (a' /a) 2 can be adjusted to op- 
timize the agreement between the physically most rele- 
vant integrals involving [/,/] and f3 JW [f, /] + ^(od v 
(V/). The simplest choice is 



(2.21) 



Of course, the IHS and EHS gases described by the 
Boltzmann equations l|2.5|l and (|2.20() . respectively, are 
intrinsically different. However, it might be possible that 
the main transport properties, which are measured by 
low- velocity moments such as in Eqs. (|2.1f) - (|2.4|l . are sim- 
ilar in both systems. As said in Sec.Q] one of the goals 
of this paper is to check this expectation in the case of 
the uniform shear flow Hal. 



III. UNIFORM SHEAR FLOW 

The uniform (or simple) shear flow (USF) is perhaps 
the nonequilibrium state most widely studied, both for 

and conventional [42j gases. In this state the gas is 
enclosed between two infinite parallel planes located at 
y = —L/2 and y = +L/2 and in relative motion with 
velocities —U/2 and +U/2, respectively, along the x- 
direction. The planes do not represent realistic bounding 
walls, in contrast to what happens in the true Couette 
flow [43. Instead, the planes represent virtual bound- 
aries where the Lees-Edwards boundary conditions are 
applied |43l 0] : every time a particle crosses one of the 
planes with a given velocity v, it is reentered at once 
through the opposite plane with a velocity v' such that 
the relative velocity with respect to the plane is pre- 
served, i.e., v' = v — Ux if the particle is reentered 
through the bottom plate and v' = v + Ux if it is reen- 
tered through the top plate. In terms of the velocity dis- 
tribution function, these generalized periodic boundary 
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conditions read 

f(y = ±L/2, v; t) = f(y = T L/2, v T U% t). (3.1) 

This process injects energy into the system. Suppose 
a particle with a velocity v crosses the top plane (i.e., 
v y > 0); it is then transferred to the bottom plane with a 
new velocity v = v — Ux. The change in kinetic energy 
is therefore proportional to v' — v 2 = 2U(U/2 — v x ), 
which is positive (negative) if v x < U/2 (v x > U/2). 
Thus, some particles gain energy while other particles 
lose energy through the boundary conditions. On the 
other hand, the shearing motion tends to produce a neg- 
ative shear stress (P xy < 0), so that particles moving up- 
ward near the top wall preferentially have v x — U/2 < 0. 
Therefore, v' — v 2 > on average. This viscous heating 
effect tends to increase the temperature, in opposition to 
the cooling effect due to the inelasticity of collisions (in 
the IHS case) or to the drag force (in the EHS case). 

Starting from a given initial condition /(r,v;0), and 
after a certain transient regime, the system is expected 
to reach a nonequilibrium steady state where the above 
heating and cooling effects cancel each other. By symme- 
try reasons, this steady state is characterized by uniform 
density and temperature, and by a linear velocity profile 
u(r) = ayx, where a — U/L represents the imposed shear 
rate. More generally, the steady distribution function be- 
comes uniform when the velocities of the particles are 
referred to the Lagrangian frame of reference co-moving 
with the flow velocity J44J: 

/(r,v)-»/(V), V = v-ayx. (3.2) 

If the initial distribution function /(r,v;0) depends on 
space only through the coordinate y normal to the plates, 
this property is maintained by the Boltzmann equations 
(1231) and (j2~2TJ|) . so that one can make 

Q 

v-V^ % -. (3.3) 

Furthermore, if /(r, v; 0) is consistent with the symmetry 
property (|3.2|l . again this is maintained by the Boltzmann 
equations (|2.5|) and (|2.2U|1 . what implies that 

W^-aF B A. (3.4) 

In this latter situation, Eqs. I|2.5[) and (|2.20l) become spa- 
tially homogeneous equations since, in agreement with 
Eq. H3.4fl . the effect of convection is played by the non- 
conservative inertial force F s h ca r = —maV y Sc. In what 
follows, we will refer to the transient solution of Eqs. 
(|2.5|) and H2.20fl with the replacement l|3.4|) as the homo- 
geneous transient problem. 

On the other hand, if the initial distribution depends 
spatially on y only but does not become uniform un- 
der the transformation l|3.2|l . then the replacement (|3.3|) 
is valid but l|3.4|) is not. In that case, Eqs. I|2.5|) and 
(|2.2U|) cannot be made equivalent to uniform equations 



and their transient solutions define the inhomogeneous 
transient problem. 

For the inhomogeneous transient problem we will con- 
sider two different initial states. The first one is the equi- 
librium state 

«"°>=^) 3/ M-^)' (3 - 5) 

where n and T° are constants (4^. The gas is initially 
at rest (in the laboratory or Eulerian frame) , but almost 
immediately the Lees-Edwards boundary conditions pro- 
duce fluid motion near the walls, this motion being sub- 
sequently propagated to the rest of the system through 
free streaming and collisions. Eventually the flow velocity 
reaches the linear profile u x (y) = ay. The transient pe- 
riod from u x — to u x — ay induces inhomogeneities in 
the density and temperature profiles, even though these 
quantities are initially uniform. As a second choice for 
the initial state we will take the distribution 

/(r, v; 0) = (|v - u°(y)| - V°) , (3.6) 

where the initial density and velocity fields are 

n%y) = n(l + ±sm^y (3.7) 

u°(y) = u(co S ^--^\% (3.8) 

respectively, while the initial temperature T° = mV° 2 /3 
is uniform. The initial state Ij3.6|) is very different from 
(|3.5|) . Now, all the particles have the same magnitude V° 
of the peculiar velocity. In addition, the initial density 
and flow velocity fields have the opposite symmetries to 
the ones imposed by the boundary conditions: n°(—y) — 
n = — \n°(y) — n] and u°(— y) = u°(y). Therefore, high 
gradients are expected during the period of time before 
the boundary conditions establish a linear velocity profile 
and uniform density and temperature. 

Regardless of the initial preparation of the system, con- 
servation of the total number of particles implies that the 
average density coincides for all times with the initial 
value n, i.e., 

1 f L/2 

n = - / dyn(y,t). (3.9) 

L J-L/2 

However, the average temperature 
1 f L ^ 2 

T(t) = — / dyn(y,t)T(y,t) (3.10) 
nL J_ L/2 

changes in time during the transient regime as a conse- 
quence of the competition between the dissipative cooling 
and the viscous heating. 
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A physically motivated way of measuring time is 
through the accumulated number of collisions per par- 
ticle s(t) from the initial state to time t. In the local 
equilibrium approximation, s(t) = So(t), where 



so(*) 



1 2 1 

2 y^n nL Jq 



J-L/2 



T(y,t>) 



(3.11) 



Here, the local characteristic time r(y, t) is given by Eqs. 
(|2.12l) and l|2.13l) . The factor | takes into account that 
two particles are involved in each collision, while the fac- 
tor 2/y / 7r accounts for the relation T m f p /r = y^r/2 [cf. 
Eq. H2.17fl ]. Note that so(t) is the (local equilibrium) 
number of collisions per particle of the IHS gas. The 
equivalent quantity for the EHS gas, s' (t), is obtained 
from Eq. (|3.11|) by using the corresponding local char- 
acteristic time T ; (y,t) defined by Eq. i|2.19|) . instead of 
r(y,t). In principle, s W ^ P a o{t) and, more gener- 
ally, s'(t) ^= /3s(t), unless the density and temperature 
profiles, and their history, are the same in both systems 
[lOj . We will come back to this point in Sec. IV Bl 

In the analysis of the homogeneous transient problem, 
we will start from a local equilibrium initial state 



/(r,v;0)=- 



3/2 



exp 



.2irT°J 

(3.12) 

In this case, the initial state is already uniform in the La- 
grangian frame [cf. Eq. 1)3.2(1 ] . so that the velocity field 
is kept linear, u x = ay, the density is constant, n = n, 
and the temperature is uniform, T(y,t) = T(t). Conse- 
quently, Eq. H3.11fl becomes 



If* dt' 
s {t) = —j= \ — — , 

V* Jo nf) 



(3.13) 



with a similar equation for s' (t) in the EHS case. Apart 
from the temperature T(t) and the elements Pij(t) of the 
pressure tensor, we will also evaluate in the homogeneous 
transient problem the ratio (V^) / (Vi 2 )o, as wen as the 
fourth and sixth cumulants, 



a 2 



(V 4 ) 



1, 



«3 



(V 6 



l + 3a 2 , (3.14) 



where (V 4 ) = 15(T/m) 2 and (^ 6 ) = 105(T/m) 3 . We 
recall that the quantity (V^ 3 2 ) is directly related to the 
cooling rate of the IHS gas via Eq. I|2.1U[) . while the 
cumulants are measures of the deviation of the energy 
distribution from the Maxwellian. Those deviations will 
also be monitored through the ratios 



Re(t) 



c: +i dw2fdvf (v,ty 



£ = 0,1,2,3. 



(3.15) 

Rt(t) is the fraction of particles that at time t move with 
a speed between We and We+i, divided by the same frac- 
tion in local equilibrium. We have taken for the integra- 
tion limits the values Wi — Cty/2T{t) Jm with Co 
d = 1, C 2 = 2, C 3 = 3, and C 4 = oo. 



0, 



In either transient problem, the final steady-state tem- 
perature T s is smaller or larger than the initial value 
T° depending on whether initially the dissipative cool- 
ing dominates or is dominated by the viscous heating, 
respectively. By dimensional analysis, T s is proportional 
to a 2 times a certain function of a, being independent 
of T°. Stated differently, the steady state is such that 

1/2 

when the steady-state characteristic time t s oc T s is 
nondimensionalized with the constant shear rate a, then 
it becomes a function of a only. More generally, the re- 
duced steady-state velocity distribution function 



1 (IT \ 3 ^ 2 V 

/ s *(c) = ±(±^) / S (V), c = ^=L=, 



(3.16) 



depends on a but is independent of the shear rate a and 
the initial preparation of the system. Owing to the sym- 
metry properties of the USF, 

fs(Cxt Cy, C Z ) — f*(C X , Cy, —C Z ) = f*( — C X , —Cy, C Z ). 

(3.17) 

Since / S *(C) depends on the three velocity components 
in a non-trivial way, it is difficult to visualize it, so that 
it is convenient to consider the following marginal distri- 
butions: 

/OO POO 
dC z dC y e(±C y )f:(C), (3.18) 
-oo J —oo 



dC z 



dC x Q(±C x )f:(C), (3.19) 



F(C) = C 2 J dC/;(C). 



(3.20) 



The function gi + \c x ) is the distribution of the x- 
component of the velocity of those particles moving up- 
ward (i.e., with C y > 0). The functions gi (C x ), 

9y + \c y ), and g y \C y ) have a similar meaning. The 
symmetry properties l(3.17|l imply that 



9y + \C y ) 



,(-) 



(3.21) 

While the functions l|3.18|l and (|3.19|l provide information 
about the anisotropy of the state, F(C) is the probability 
distribution function of the magnitude of the peculiar 
velocity (in units of the thermal speed), regardless of its 
orientation. 



IV. MONTE CARLO SIMULATIONS 

We have solved numerically the Boltzmann equation 
(|2.5|) for the IHS system and the Boltzmann equation 
(|2.20ll for the EHS system by means of the Direct Simu- 
lation Monte Carlo (DSMC) method |H El E3 . For the 
sake of completeness, we give below some details about 
the application of this method to our problem. 
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A. Inhomogeneous transient problem 

The USF has been implemented in the inhomogeneous 
transient problem by applying the Lees-Edwards bound- 
ary conditions (|3.1|) , using the form l|3.3[) for the convec- 
tion operator, and starting from the initial distribution 
()3.5|1 or (|3.6|1 . The separation between the boundaries 
has been taken as L = 2.5A and the shear rate has been 
fixed at ar° = 4, where 



A 



(4.1) 



are the initial (global) mean free path and characteris- 
tic time, respectively, of the IHS gas. The coefficient of 
restitution for the IHS has been taken as a = 0.9. This 
same value has been taken in the EHS case for the fric- 
tion constant 7 = ^Co and the factor /3, as given by Eqs. 
(I2.16|) and l|2.21|l . respectively. 

According to the DSMC method HEH^, the system is 
split into M layers of width Sy = L/M. The velocity dis- 
tribution function is represented by the positions {yi(t)} 
and velocities {v^(i)} of a set of TV simulated particles: 



/(l/,v;t) 



1 N 



W(t))«(v-Vi(i)), (4.2) 



where A = (N/L) jn is a constant formally representing 
the area of a section of the system normal to the y-axis, 
so that ASy represents the volume of a layer. The number 
of particles inside a given layer I = 1 , . . . , M is 



JV 



(4.3) 



where &i(y) is the characteristic function of layer I, i.e., 
&i(y) = 1 if y belongs in / and 0/(y) = otherwise. The 
(coarse-grained) number density, mean velocity, temper- 
ature, and pressure tensor of layer / are 



m(t) = 



Ni(t) 
ASy ' 



(4.4) 



Ul(t) = jvTi)? e/(yi( * ))v<w ' (4 - 5) 



N 



Tl(t) = m(t) E edm{t)) [Mt) Ul(t)]2 ' (4 ' 6) 



JV 

p ' w = E e ^)) [ v * (*) - u ^wi ^(i) - u z (t)] . 

(4.7) 



The average temperature and pressure tensor along the 
system are given by 



1 M 



1=1 



1 M 



(4.8) 



(4.9) 



1=1 



The positions and velocities {vi(t)} of the par- 

ticles are updated from time t to time t + St in two stages: 



1. Free streaming. In this stage, 

y l {t + 8t) = yi(t)+v iy (t)St. 



(4.10) 



If particle i crosses the top wall, i.e., yi(t + 5t) > 
L/2, then its position and velocity are redefined as 



yi{t + 6t)->yi(t + 6t)-L, Vi(t) -»• Vi(t)-aL5L (4.11) 

A similar action takes place if particle i crosses the 
bottom wall, i.e., yi(t + St) < —L/2. 

In the case of IHS, the velocities are not modi- 
fied during the free streaming stage. In the case 
of EHS, however, the action of the (local) friction 
force yields 

Vi(t + St) = u/(t) + [v<(t) - uj(t)] e^' WSt/2 . (4.12) 

Here, / is the layer where particle i sits at time 
t and Ci(t) n i{t)\jTi{t)(l — a 2 ) is the coarse- 
grained version of the cooling rate Co defined by 
Eq. (|2~To|) . 



2. Collision stage. In this stage, a number 



2V2N/M A° 



(IHS), 



A/} - /3(a) 



wjSt 



2V2N/M A 



(EHS) 



(4.13a) 



(4.13b) 



of candidate pairs are randomly selected for each 
layer /. In Eqs. (|4.13|) . A is given by Eq. I|4.1|) 
and u>/ oc ^jTjJt) is an upper estimate of the max- 
imum relative speed in layer /. The collision be- 
tween each candidate pair ij is accepted with a 
probability equal to Vij/wi, where Vij is the rel- 
ative speed. If the collision is accepted, a direction 
<r is chosen at random with equiprobability and the 
velocities (vj,v_y) are replaced by (v£,v^), accord- 
ing to the collision rule (|2.8|l with a < 1 (IHS) or 
a = 1 (EHS). 
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The numerical values for the "technical" parameters 
are as follows. The layer thickness is Sy — 0.05 A (i.e., 
the number of layers is M = 50), the time step is 

St = 10- 3 t°^Jt°/T (so that it changes in time as the 

global characteristic time r oc 1/Vt does), and the to- 
tal number of particles is N = 10 4 . The hydrodynamic 
quantities (density, flow velocity, temperature, and pres- 
sure tensor) are updated every 4 time steps and recorded 
every 160 time steps. Moreover, in order to improve 
the statistics, all the quantities are further averaged over 
three independent realizations of the system. 



B. Homogeneous transient problem 

In the homogeneous transient problem the USF is im- 
plemented by working directly in the Lagrangian frame of 
reference [cf. Eq. (|3.2[l ] and using l|3.4fl . Since the result- 
ing Boltzmann equation is uniform, only the (peculiar) 
velocities {Vi(i)} of the N particles need to be stored 
and there is no need of splitting the system into cells or 
applying boundary conditions. The velocity distribution 
function is described by 



1 N 



(4.14) 



where the constant fl = N/n formally represents the vol- 
ume of the system. The temperature and the pressure 
tensor are evaluated as 



N 



P W = S5> i(<)VjW ' 



N 



(4.15) 



(4.16) 



Analogously, the fourth and sixth cumulants of the dis- 
tribution function [cf. Eq. I|3.14[l ] and the ratios Ri [cf. 
Eq. 1|3.15[) ] are computed as 



a 2 (t) 



m 1 



15T 2 (i) ^ 



E^)- 1 . ( 4 - 17 ) 



03 (t) = - 



777 1 



N 



105T3(i) TV ^ 



]Tvf(i) + l + 3 a2 (t), (4.18) 



of the relative speed, a sample of N p pairs is randomly 
chosen out of the total number N(N — l)/2 of pairs, so 
that 



(n 3 2 >w 



-E 



|V 2 (t)-V,(t)| 3 . 



(4.20) 



The velocity update {Vi(t)} — > {Vi(t+5t)} takes place 
again in two stages. In the free streaming stage for IHS, 
only the x-component of the velocities change according 



to the inertial force F s i 1C ar = — maV y x: 



V^t + 5t) = V lx {t) - V iy (t)aSt. 



(4.21) 



On the other hand, for EHS we have the drag force 
Fdrag = — m(£b/2)V in addition to F s h oa r, so that 



V ix (t + St) = [V lx {t) - V ly {t)aSt] e-^') 5 '/ 2 , 



(4.22) 



V ly At + St) = V W:Z (t)e- io{t)5t/2 , 

where Co(i) is given by Eq. (|2.16l) with r(t) = 
\°/^/2T(t)/m. 

The collision stage proceeds essentially as in the in- 
homogeneous case, except that formally the number of 
layers is M = 1. Therefore, a number 



Af=A=^ (IHS), 



AT = /3(a) 



2V2 A° 



N wSt 



(EHS) 



(4.23a) 



(4.23b) 



of candidate pairs are randomly select ed ou t of the total 
number of pairs in the system, w oc \/T{t) being an up- 
per estimate of the maximum value of the relative speeds 
{Vij} in the whole system. 

In the simulations of the homogeneous transient prob- 
lem we have considered two values of the shear rate 
(ar° = 4 and ar° = 0.1) and ten values of the coefficient 
of restitution (a = 0.5-0.95 with a step Aa = 0.05), 
both for IHS and EHS. This gives a total of forty dif- 
ferent systems simulated. However, in the discussion of 
the transient problem we will mainly report results cor- 
responding to a = 0.5 and a = 0.9. Once a steady state 
is reached, its properties are obtained by averaging the 
fluctuating simulation values over time. In the analysis 
of the steady state all the values a = 0.5-0.95 will be 
considered. 

The technical parameters of the homogeneous simu- 
lations arc St = 10~ 3 r° y^T° /T (the average quantities 
being updated every 4 time steps and recorded every 160 
time steps), N = 10 4 , and N p = 2.5 x 10 4 . 



Ri(t) 



N- 1 



N 



V(C e+1 ) - V(C e ) 

xe(w i+1 (t)-Vi{t)), 



5>(t$(t)-wi(t)) 



(4.19) 



where V(x) = erf(cc) — 2xe x / t/tt, erf (a;) being the error 
function. In order to evaluate the average of the cube 



V. RESULTS 

In this Section we present the simulation results ob- 
tained for the main physical quantities in the USF and 
compare the properties of the genuine IHS gas with those 
of the EHS gas. 



8 




FIG. 1: (Color online) Accumulated number of collisions per 
particle as a function of time for IHS [s(t), smooth solid line] 
and EHS [s'(t), smooth dashed line] in the case a = 0.9, ar° = 
4. The fluctuating solid line represents the ratio s' (t) / s{t), the 
corresponding scale being that of the right vertical axis. The 
data have been obtained starting from the initial distribution 
function 13.51 1. 



A. Inhomogeneous transient problem 

As said in Sec. IIVI the initial distribution function in 
the inhomogeneous transient problem is either that of 
equilibrium, Eq. I|3.5[) . or a strongly nonequilibrium one, 
Eq. (|3.6[) . We have restricted ourselves to a coefficient 
of restitution a = 0.9 and a shear rate a = 4/r° = 
(4/0.95)/t ', where r° and r ' = r°//3(a) = r°/0.95 are 
the initial (global) characteristic times of the IHS and 
EHS gases, respectively. The values of a and a are such 
that the viscous heating initially prevails over the dissi- 
pative cooling (either collisional or frictional) and so the 
temperature increases and the mean free time decreases. 

We will mainly monitor the temporal evolution of the 
physical quantities by using an internal clock, namely the 
accumulated number of collisions per particle s(t) (IHS) 
and s'(t) (EHS), rather than the external time t. The 
quantities s(t) and s'(t) are computed directly by divid- 
ing the total number of accepted collisions until time t 
by the total number of particles; in general, they slightly 
differ from the local equilibrium values so(t) and s' (t) 
[cf. Eq. (|3.11|) ]. Insofar as the velocity distribution func- 
tion /(y,v;i) is similar in both systems, one can expect 
that s'(t)/f3(a) ~ s(t) jib]. Figure d shows s(t) and s'(t) 
as functions of time in the case of the initial state l|3.5[l . 
The slopes of those curves are directly related to the re- 
spective temperatures; both slopes increase monotoni- 
cally until becoming constant for t/r° > 8. As expected, 
the accumulated number of collisions in the EHS gas up 
to any given time t is smaller than in the IHS gas, i.e., 
s'(t) < s(t). Figure ^ also shows the temporal evolution 
of the ratio s'(t)/s(t). It fluctuates around (3(a) = 0.95 
up to t/r ~ 4 and stays very close to 0.95 thereafter. 
This already provides an indirect validation of the prac- 
tical equivalence between the profiles and history of the 
hydrodynamic quantities in both systems. 

Figure [21 shows the velocity, density, and temperature 



profiles at times t/r° — 0.13, 0.5, 1, 1.5, and 2.0. The cor- 
responding accumulated numbers of collisions per parti- 
cle (for IHS) are s ~ 0.08, 0.4, 1.1, 2.1, and 3.4 in the case 
of the initial distribution (|3.5|) . and s ~ 0.09, 0.5, 1.5, 2.8, 
and 4.6 in the case of the initial distribution (|3.6() . By 
time t/r — 0.13 only about 16-18% of particles have 
collided, so that the deviations from the initial profiles 
are essentially due to the boundary conditions. As a 
consequence, at t/r° = 0.13 the flow velocity is almost 
the initial one everywhere except in the layers adjacent 
to the walls; moreover, those layers have a much larger 
temperature than the bulk, while the number density is 
still practically unchanged. As time advances, the more 
energetic particles near the boundaries travel inside the 
system and transfer part of their momentum and energy 
to the other particles by means of collisions. This pro- 
duces a stretching of the shape of the velocity profile as 
well as a homogenization of temperature. The layers ad- 
jacent to the walls are first partially depopulated in favor 
of the central layers, but as the velocity profile becomes 
linear and the temperature becomes uniform, so does the 
density. In summary, Fig.[2clearly shows that the hydro- 
dynamic profiles freely evolve toward the characteristic 
profiles of the USF, regardless of the initial preparation 
of the system. It must ne noted that, although the USF 
is known to be unstable with respect to excitations of suf- 
ficiently long wavelengths [T^. l2ll I23L |29| , the size of the 
simulated systems (L — 2.5A ) is small enough to sup- 
press such an instability. As a matter of fact, a recent 
analysis from kinetic theory shows that at a = 0.9 
the instability does not appear unless L > 25 A . 

For the sake of clarity of the graphs, in the remainder 
of this Subsection we restrict ourselves to present simu- 
lation results obtained from the initial distribution <|3.5[) . 
To monitor the time needed to establish a linear veloc- 
ity profile and uniform density and temperature, we have 
followed the evolution of those quantities averaged over 
the four central layers (— 2Sy < y < 2Sy) and over the 
three top layers (L/2 — 3Sy < y < L/2). The results are 
displayed in Fig. |31 The flow velocity at y fts fluctuates 
around zero, as expected by symmetry, whereas the ve- 
locity near the top wall monotonically increases (except 
for fluctuations) toward the wall velocity U/2. The max- 
imum density difference appears after about one collision 
per particle, while the maximum temperature difference 
appears earlier at s ~ 0.15. The characteristic hydrody- 
namic profiles of USF, i.e., linear velocity and uniform 
density and temperature, are reached approximately at 
s = 4 (which corresponds to a "real" time t/r ~ 2.2). 
From this time on, the evolution proceeds essentially as 
in the homogeneous transient problem. 

In Figs. [21 and [3] we have scaled the local tempera- 
ture T(y,t) with respect to its average value T(t) in 
order to focus on the transient period toward unifor- 
mity. However, once the system becomes "homoge- 
neous" [in the sense of Eq. |3.2p ] at s ~ 4, the tem- 
perature keeps evolving in time until the steady state 
is reached. This is observed in the top panel of Fig. 
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FIG. 2: (Color online) Hydro dynamic profiles for IHS (solid lines) and EHS (dashed lines) in the case a = 0.9, ar = 4, at 
times t/r° = 0.13, 0.5, 1, 1.5, and 2. The data have been obtained starting from the initial distribution function 13.51 (curves 
without symbols) or 13.01 (curves with symbols). 



which shows the evolution of T/T°. We note that 
the total transient period is much longer than the du- 
ration of the inhomogeneous state. The temperature 
reaches a stationary value much larger than the initial 
one (T s /T° ~ 147) at s ~ 50 (t/r ~ 9.8). Station- 
arity of temperature does not necessarily imply that the 
steady state has been reached since, in principle, the par- 
ticles could redistribute their velocities along time with- 
out altering the mean kinetic energy. A strong indication 
that this is not actually the case is provided by the bot- 
tom panel of Fig. which shows the evolution of the 
(reduced) shear stress —P xy /nT and the (reduced) nor- 
mal stress difference (P xx — P yy )/nT. These two quan- 
tities reach stationary values —P X y,s/nT s ~ 0.33 and 



(Pxx,b - P yy ,s)/nT s ~ 0.24 after s ~ 40 (t/r ~ 8.3). 

In the above comments on Figs.|2JE]we have focused on 
the physical features of the transient toward the steady- 
state USF, without distinguishing between the IHS and 
EHS systems. In fact, as Figs. show, the results con- 
cerning the hydrodynamic quantities and their fluxes are 
practically identical in both systems, even when high gra- 
dients are transitorily present. Therefore, the transport 
properties of a gas of IHS can be satisfactorily mimicked 
by a gas of EHS having an adequate diameter and sub- 
ject to an adequate friction force. However, although a 
coefficient of restitution a = 0.9 is rather realistic, it is 
reasonable to expect that this approximate equivalence 
IHS^EHS deteriorates as dissipation increases. To ana- 
lyze this expectation we have considered other values of 
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FIG. 4:JColor _online)_ Evolution of T/T° (top panel) and 
-Pxy/nT and (P xx - P yy )/nT (bottom panel) for IHS (solid 
lines) and EHS (dashed lines) in the case a = 0.9, ar° = 4. 
Time is measured by the accumulated number of collisions 
per particle (s) in the case of IHS and by the same quantity, 
but divided by (5, (s'/f3) in the case of EHS. The data have 
been obtained starting from the initial distribution function 
11531 . 



FIG. 3: (Color online) Evolution of the flow velocity, the den- 
sity, and the temperature around y = and y — L/2 — (3/2)Sy 
for IHS (solid lines) and EHS (dashed lines) in the case 
a — 0.9, ar° = 4. Time is measured by the accumulated 
number of collisions per particle (s) in the case of IHS and 
by the same quantity, but divided by /3, (s'//3) in the case of 
EHS. The data have been obtained starting from the initial 
distribution function 13.51 . 

a in the homogeneous transient problem and, especially, 
in the steady-state properties. 

B. Homogeneous transient problem 

As discussed in Sec. IIIII the Boltzmann equation for 
USF allows for solutions which are spatially uniform 
when the velocities are referred to the (local) Lagrangian 
frame of reference [cf. Eq. I|3.2|) ], The simulation of 
the corresponding Boltzmann equation by the DSMC 
method is much simpler than in the inhomogeneous prob- 
lem, as described in Sec. IIVBI In these homogeneous 
simulations we have considered a = 0.5-0.95 (with a step 
Aa = 0.05) and two values of the shear rate, namely 
cit° = 4 and ar° = 0.1. The former value is large enough 
to make viscous heating initially dominate over (inelastic 
or frictional) cooling, even for a = 0.5, so that T(t) > T°. 



Conversely, ar° — 0.1 is small enough to produce the op- 
posite effect, T(t) < T°, even for a = 0.95. On the 
other hand, at a given value of a, the intrinsic steady- 
state properties must be independent of the value of the 
shear rate, and this will provide an important indicator 
to determine whether the steady state has been reached 
or not. 

Although we have performed simulations for the ten 
values a = 0.5-0.95, in this Subsection we will focus on 
two values: a = 0.9 (moderately small dissipation) and 
a = 0.5 (large dissipation). In addition to the simula- 
tion data for IHS and EHS, we will present the results 
from the solution of an extension 0,0] of the Bhatnagar- 
Gross-Krook (BGK) model which is inspired in the 
approximate equivalence IHS<->EHS. The solution of this 
BGK-like model for the USF problem is worked out in 
Refs. HSIH. 

We begin by showing the accumulated number of col- 
lisions per particle as a function of time in Fig. where 
the case a = 0.7 is also included. As said in connection 
with Fig.Q] the slope of s(t) and s'(t) is proportional to 
y/T(t). At a shear rate ar° = 4, viscous heating dom- 
inates and so the temperature monotonically increases, 
especially for a — 0.9; on the other hand, at ar° = 0.1, 
dissipative cooling prevails and so the temperature mono- 
tonically decreases, especially for a — 0.5. The almost 
perfect agreement between s(t) and s'(t)/j3 for a = 0.9 
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FIG. 5: (Color online) Accumulated number of collisions per 
particle as a function of time for IHS [s(t), solid lines] and 
for EHS [s'(t), dashed lines], in the latter case divided by 
/3 = |(1 + a), in the cases a = 0.5, 0.7, and 0.9. The left 
panel corresponds to a shear rate ar° = 4, while the right 
panel corresponds to ar° = 0.1. The dotted lines are the 
predictions obtained from the solution of a BGK model. Note 
that in the right panel the curves corresponding to IHS, EHS, 
and BGK at a — 0.9 are practically indistinguishable. The 
data have been obtained starting from the initial distribution 
function H3.12H . 



is an indirect indication that T(t) is practically the same 
for the IHS gas and the EHS gas, as will be confirmed 
later on. However, as the inelasticity increases, so does 
the temperature difference in both systems, the EHS sys- 
tem having a slightly higher temperature than the IHS 
system at any given time t. Despite its simplicity, the 
BGK model does quite good a job, but it tends to over- 
estimate both s{t) and s'(t)/(3. It must be said that the 
BGK curves for s(t) have actually been obtained from the 
local equilibrium approximation i|3.13|) . As we will see, 
the BGK temperature presents a very good agreement 
with the simulation results for EHS, so that the small 
discrepancies between the DSMC curves for s'it)/ (i and 
the BGK curves are essentially due to the local equilib- 
rium approximation s(t) — ► srj(t) used in the latter. 

To confirm this point, we plot in Fig. EJthe ratio be- 
tween the actual number of accumulated collisions per 
particle and the local equilibrium estimate obtained from 
Eq. I|3.13[) by a numerical integration using the actual 
values of temperature. Except for a short initial pe- 
riod, s(t)/so(t) and s'(t)/s' Q (t) take values smaller than 
1 and tend to steady-state values practically indepen- 
dent of the shear rate. The instantaneous collision rate 
is proportional to the average value of the relative speed 
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FIG. 6: (Color online) Evolution of the ratio between the 
actual number of collisions per particle and the local equilib- 
rium estimate for IHS (solid lines) and EHS (dashed lines) in 
the cases a = 0.5,0.9 and ar° = 0.1,4. The data have been 
obtained starting from the initial distribution function 13.121 . 



(V12), which in the local equilibrium approximation is 
(V12) -» (^12)0 = y/2(V) = i^T/i^. Thus, the 
fact that s(t) < so(t) is consistent with a noncquilib- 
rium velocity distribution such that the mean speed is 
smaller than the local equilibrium one, i.e., (V) < (V)o, 
this effect being more noticeable as the dissipation in- 
creases. By definition, (V 2 ) = (V 2 )o = 3T/m. There- 
fore, the inequality (V) < (V)q indicates an underpop- 
ulation of the region of moderately low velocities of the 
nonequilibrium distribution function (with respect to the 
Maxwcllian) that compensates for an overpopulation of 
the high-velocity region, as will be seen later. 

The evolution of the relative temperature T/T°, the 
reduced shear stress —P xy /nT, and the reduced normal 
stress difference {P xx — P vy )/nT is displayed in Figs. 
and [3] for a = 0.9 and a = 0.5, respectively. In the for- 
mer case, an excellent agreement between the simulation 
results for both types of system exists. In addition, the 
theoretical results obtained from the BGK model accu- 
rately describe the behavior of the simulation data. In 
the case a = 0.5, however, the EHS system tends to have 
a larger temperature than the IHS system, the steady- 
state value being about 12% larger in the former case 
than in the latter. This is partially due to the fact that 
the true cooling rate £ of the IHS gas [cf. Eq. (|2.10Jl ] is 
larger than the local equilibrium value Co [cf- Eq. 12.16fl ] 
imposed on the EHS gas, as we will see later on. This 
also explains why the BGK model, which also makes use 
of the approximation £ — * Co , predicts a temperature in 
good agreement with the simulation data for EHS. Since 
the EHS temperature is larger than the IHS one, the 
shear rate normalized with the collision rate is smaller 
for EHS than for IHS and, consequently, the distortion 
with respect to the Maxwellian, as measured by the shear 
stress and, especially, by the normal stress difference, is 
smaller in the former case than in the latter. Compar- 
ison between Figs. and |H1 shows that the duration of 
the transient period (as measured by the number of col- 
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FIG. 7: (Color online) Evolution of T/T°, -P xy /nT, and 
(Pxx - P vy )/nT for IHS (solid lines) and EHS (dashed lines) 
in the cases a = 0.9 with ar° = 0.1 and ar° — 4. The dotted 
lines are the predictions obtained from the solution of a BGK 
model. Note that in the top panel the curves corresponding 
to IHS, EHS, and BGK are practically indistinguishable. The 
data have been obtained starting from the initial distribution 
function 13.121 1. 



lisions per particle) decreases as the inelasticity increases: 
s ~ 40 at a = 0.9 versus s ~ 20 at a = 0.5. However, 
when that duration is measured in real units, it depends 
mainly on the shear rate and not on a, namely t ~ 9r° 
for ar° = 4 and t ~ 160-200r° for ar° = 0.1. 

It is interesting to note the similarity between the 
curves in Fig. 0] and those corresponding to ar° = 4 in 
Fig. This is made clear in Fig. |5J where the evolution 
of the global quantities in the IHS case for the homoge- 
neous and inhomogeneous transient problems are shown. 
In the three situations considered, the initial shear stress 
and normal stress differences are zero; however, given 
the high value of the shear rate, the velocity distribution 
function rapidly changes and adapts itself to the imposed 
shear rate, giving rise to a sharp maximum of the reduced 
shear stress and normal stress difference. Henceforth, 
as the temperature increases so does the collision rate, 
so that the relative strength of the shear rate becomes 
smaller and smaller until the steady state is reached at 




FIG. 8: (Color online) Evolution of T/T°, -P xy /nT, and 
(Pxx - P yv )/nT for IHS (solid lines) and EHS (dashed lines) 
in the cases a — 0.5 with ar° — 0.1 and ar° — 4. The dotted 
lines are the predictions obtained from the solution of a BGK 
model. The data have been obtained starting from the initial 
distribution function 13.12H . 



s « 40. The first stage lasts about one collision per par- 
ticle, has a kinetic nature, and is sensitive to the initial 
preparation of the system. On the other hand, the subse- 
quent relaxation toward the steady state defines a much 
longer hydrodynamic stage that becomes more and more 
independent of the initial state, provided that T° < T s . 
A similar hydrodynamic regime exists for the class of ini- 
tial states with T° > T s (e.g., for ar° = 0.1). These 
comments also apply to the two cases with a = 0.5. 

In general, at a given value of a, the intrinsic veloc- 
ity distribution function in the hydrodynamic regime de- 
pends on time only through its dependence on the shear 
rate nondimensionalized with the (time-dependent) col- 
lision rate. More specifically, 



/(V,<) = 



2T(t) 



3/2 



f*(C(t);a*(t)), (5.1) 



where 



C(t) = 



V 



y/2T(t)/m 



<(t)=aT v (t). (5.2) 
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FIG. 9: J Color online) JEvolution of T/T° (top panel) and 
—Fxy/r^T and (P xx — Pyy)/nT (bottom panel) for IHS in 
the case a — 0.9, ar° — 4. The dashed and dash-dot lines 
correspond to the inhomogeneous transient problem with the 
initial conditions (13.511 and 13, 6L respectively, while the solid 
lines correspond to the homogeneous transient problem with 
the initial conditions 13,1211 . 



The definition of the reduced shear rate a* by means of 
the viscosity characteristic time t 1] instead of the mean 
free time r m ft or the characteristic time r is an irrelevant 
matter of choice since these three quantities differ only 
by constant factors, namely T m f t /r^ ~ 0.787 and t/t 71 ~ 
0.888. In the steady state, a*\t ) -> a* and f*(C;a*) -> 
/*(C;<)=/*(C) [cf.Eq.E3SJ]- 

During the hydrodynamic relaxation toward the steady 
state, it is insightful to define a time-dependent shear 
viscosity r](t) = —P xy (t)/a. As a consequence of Eq. 
(|5.1|) . the ratio r)(t)/r)o(t), where f]o(t) = nT(t)T v (t) is 
the Navier-Stokes viscosity in the elastic limit [cf. Eq. 
ipngfr]. depends on time only through a nonlinear de- 
pendence on a*(t). Figure ITUI shows the reduced shear 
viscosity 77* = 77/770 versus the reduced shear rate a* in 
the case a — 0.5. The values plotted correspond to a 
temporal window 2 < s < 50. The lower limit guarantees 
that the system has practically lost memory of its initial 
state, while the upper limit is long enough to guarantee 
that the steady state (represented by an open symbol) 
has been reached. For each case (IHS, EHS, or BGK), 
the steady-state point (a*, 77* ) splits the respective curve 
into two branches: the one to the left of the point corre- 
sponds to T < T° (e.g., ar° — 0.1), while the branch to 
the right corresponds to T > T° (e.g., ar° = 4). The en- 
tire curve represents the nonlinear shear viscosity 77* (a*) 
for a = 0.5, and the steady-state point 77* = 77* (a*) is just 
a particular (and singular) value [6J. The formal extrap- 
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FIG. 10: (Color online) (Transient) reduced shear viscosity 
rj*(t) = T](t)/r)o(t) versus the (transient) reduced shear rate 
a*(t) = aT v (t) for IHS (solid lines) and EHS (dashed lines) 
in the case a = 0.5. The dotted lines are the predictions ob- 
tained from the solution of a BGK model. The circle, square, 
and diamond are the steady-state points for IHS, EHS, and 
BGK, respectively. In each case, the curve to the left of the 
symbol corresponds to ar° = 0.1, while the curve to the right 
of the symbol corresponds to ar° = 4. 



olation of 77* (a*) to zero shear rate gives the (reduced) 
Navier-Stokes shear viscosity 77^0 = lim Q *_»o 7 7*( a *) a t 
a = 0.5. The expected values are "0,01 77^,3 — 1.3 for IHS 
and 77^3 ~ 1.1 for EHS and BGK. It is worthwhile noting 
that, except for fluctuations in the simulation data, the 
curves 77* (a*) for IHS, EHS, and BGK practically coin- 
cide, at least in the interval 0.2 < a* < 1.6. Thus, the 
main difference among the three approaches lies in the 
steady-state point where the system "decides" to stop. 
For a more extensive discussion on the rheological func- 
tion 77* (a*) and the distinction between 77* and 77^ the 
reader is referred to Ref. 

Although the most relevant information in the USF 
problem is conveyed by the elements of the pressure ten- 
sor, they of course do not exhaust the physical informa- 
tion one can extract from the simulations. In Figs. 1111 
and 1 121 we present the evolution of the ratio (V^)/ (^12)0 
(where V12 is the relative speed between a pair of parti- 
cles) and the fourth and sixth cumulants [cf. Eq. I|3.14|l ] 
for a = 0.9 and a = 0.5, respectively. The average value 
(V12) is physically interesting because it is proportional 
to the true cooling rate ( of the IHS gas. In fact, we have 
checked in the simulations that the value of £ obtained 
from Eq. p.lOfl agrees with the one computed directly 
from the collisional energy loss, Eq. I|2.9|l . However, the 
second method is much noisier than the first one since 
it involves only the colliding pairs, whereas all the pairs 
contribute to (V^) [see, however, Eq. (|4.20|) and the com- 
ment above it]. 

As happened with Figs.|7|and|S| the simulation data for 
each one of the three plotted quantities in the case ar° = 
4 present a respective maximum during the kinetic stage 
of the evolution, representing the largest departure from 
the Maxwellian. It is noteworthy that the fluctuations of 
the quantities plotted in Figs.llllandll2lare much larger 
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FIG. 11: (Color online) Evolution of the ratio (V^) / (V^o, 
the fourth cumulant 0,2, and the sixth cumulant —az for IHS 
(solid lines) and EHS (dashed lines) in the cases a — 0.9 with 
ar° = 0.1 and ar° — 4. The data have been obtained starting 
from the initial distribution function 13.12H . 
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FIG. 12: (Color online) Evolution of the ratio (V? 2 ) / {V? 2 )o, 
the fourth cumulant a 2 , and the sixth cumulant —03 for IHS 
(solid lines) and EHS (dashed lines) in the cases a = 0.5 with 
ar° = 0.1 and ar° = 4. The data have been obtained starting 
from the initial distribution function 13.121 . 



in the IHS case than in the EHS case. Otherwise, the 
temporal evolution and the steady-state values are very 
similar in both systems. 

We observe that (V^) > (^12)0, the relative differ- 
ence increasing with the inelasticity. This explains that 
the (internal) collisional cooling rate £ of the IHS gas is 
larger than the (external) frictional cooling rate £0 im- 
posed on the EHS gas, as said before in connection with 
Fig. |H1 However, even though the imposed cooling rate 
for EHS is the local equilibrium one, this system satis- 
factorily mimics the distortion from local equilibrium, as 
measured by (V± 2 } / (Vi 2 )o, of the IHS system. 

The cumulants basically probe the high-energy tail of 
the distribution, especially in the case of 03. The posi- 
tive value of the fourth cumulant 02, both for IHS and 
EHS, is a reflection of a strong high-energy overpopu- 
lation (with respect to the Maxwellian) induced by the 
shearing. This overpopulation effect is larger than and 
essentially independent of the one typically present in 
homogeneous states of granular gases jUil, which is 
absent in the EHS gas. The steady-state value of the 




5 10 15 20 

FIG. 13: (Color online) Evolution of the ratios Ro, Ri, R2, 
and R 3 [cf. Eq. lETTSl l for IHS (solid lines) and EHS (dashed 
lines) in the case a — 0.5 with ar° = 4. The data have been 
obtained starting from the initial distribution function 13.121 . 



15 



sixth cumulant 03 is practically zero for a = 0.9. On the 
other hand, for a = 0.5 one has — a 3 > 0, what is again 
related to the high-energy overpopulation. 

A more direct information about the evolution of the 
velocity distribution is provided by the ratios Ri defined 
by Eq. I|3.15() . They are plotted in Fig. ^] for the case 
a = 0.5 and ar° = 4. The maximum deviation from 
the Maxwellian takes place during the kinetic stage (s < 
1). Thereafter, the curves smoothly relax toward their 
steady-state values R ~ 1.18, R\ ~ 0.81, R2 — 1.39, 
and i?3 ~ 15, the relaxation times being practically the 
same (s w 10) in the four cases. During the hydrody- 
namic transient regime and in the steady state, the pop- 
ulation of particles moving with a speed larger than three 
times the thermal speed Vo(t) = ^j2T(t)/m is remarkably 
larger than the one expected from a Maxwellian distri- 
bution. This overpopulation effect induced by shearing 
is also present in the interval 2vo(t) < V < 3v n (t). In the 
hydrodynamic transient regime, between 90% and 93% 
of the particles move with a speed smaller than 2vo(t) 
(in contrast to the equilibrium value of 95.4%) and this 
is then the relevant region for the low-degree moments. 
While the low- velocity region V < Vo(t) is about 20% 
overpopulated with respect to the Maxwellian, the inter- 
mediate region v (t) < V < 2vo(t) is underpopulated by 
about the same amount. In fact, the fraction of particles 
moving with V < v (t) and v (t) < V < 2v (t) is about 
51% and 42%, respectively, in the steady state, while the 
corresponding equilibrium values are 42.8% and 52.6%, 
respectively. As Fig.ll3lshows. all these features are suc- 
cessfully captured by the EHS system, even at this rather 
high inelasticity. 

The fact that (V^) > (Vi 2 )o suggests that a better 
agreement between the dynamics of the IHS gas and that 
of the EHS gas could be expected if the friction con- 
stant of the latter were not chosen as 7 = |£q but as 
7 = ICoO'il) / (^2)0, so that the cooling rate of the EHS 
gas would be exactly the same functional of / as that of 
the IHS gas. To test this expectation, we have carried 
out supplementary simulations of the EHS system with 
this more refined friction constant in the case a = 0.5, 
ar° = 4. The corresponding curves are not included in 
Figs. El El H CHI 112 and US for the sake of clarity. The 
results show that, in general, the quantities associated 
with the low-order moments (temperature and pressure 
tensor) are indeed closer to the IHS values than with the 
simpler choice 7 = |Co- As a matter of fact, the EHS 
temperature is now slightly smaller, instead of slightly 
larger, than the IHS temperature. However, in the cases 
°f (^L2)/(^L2)o, 02, and 03, the results obtained in the 
EHS simulations with 7 = ^Co^ilV^i^o turn out to be 
not necessarily better than those obtained with 7 = |£o- 



C. Steady state 

Once we have analyzed the transient period toward the 
steady state for the representative cases a = 0.5 and a = 



1.0 
0.8 
0.6 
* 0.4 
0.2 

0.5 
0.4 

s"o.3 



0.2 
0.1 



0.8 
0.6 



0.2 



0.06 



X 0.04 

a* 

' P 0.02 



0.00 



-* * * 

*• 3.... 

O IHS, a=0.1/x° 
V IHS, a=4/T° 
□ EHS, a=0.1/x° 
A EHS, a=4h° 
-■■ BGK , 



•8 8 



8 



•fi- 



•ft.. 



■tS:. 9 



a. 



0.5 0.6 0.7 OS 

a 



0.9 1.0 



FIG. 14: (Color online) Steady-state values of the reduced 
shear rate a* = ar v , the reduced shear stress —P xy /nT, and 
the reduced normal stress differences (P xx — P yy )/nT and 
(Pzz — P yy )/nT as functions of the coefficient of restitution a. 
The open symbols are simulation results for IHS and EHS and 
two values of the shear rate, while the dotted lines correspond 
to the solution of a BGK model. Note that in the latter model 
Pzz,s = Pyy,s- The filled diamonds represent simulation data 
of EHS with a friction constant 7 = | (0 {V12} / (V12 )o in the 
case a = 0.5, ar° = 4. 



0.9, let us report on the most relevant steady-state prop- 
erties for all the values a = 0.5-0.95 we have considered. 
The quantities associated with the second-degree velocity 
moments, namely the reduced shear rate a* = aT n (T s ) oc 
l/v 7 ?!, the reduced shear stress —P XVtS /nT s , and the re- 



duced normal stress differences {P x 



P yy ,s)/nT s and 



(Pzz s — Pyy,s)/ n T s , are plotted in Fig. The over- 
lapping of the data obtained from simulations with the 
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two different values of the shear rate (ar° = 0.1 and 
(it q = 4) confirms that the steady state has actually been 
reached and that the intrinsic velocity distribution func- 
tion /*(C) [cf. Eq. (|3.16|) ] depends only on a and not on 
the initial preparation of the system. For a > 0.7 there 
exists a good agreement between the EHS and IHS re- 
sults for the quantities a* and —P xy , s /nT s , which are the 
most relevant properties in the USF problem. For larger 
inelasticity, however, the steady-state temperature T s is 
larger for the EHS gas than for the IHS gas, and so a* 
is smaller in the former case than in the latter. This im- 
plies that the departure from isotropy is slightly smaller 
in the EHS gas than in the IHS gas, and so is the shear 



stress 

(P X X,i 



-Pxy,s/nT s . As for the normal stress differences 



Pyy,s)/nT s and, especially, (P z 



Pyy,s)/ n Ts, 



they start to differ in both systems for a < 0.85. It is 
worthwhile noting that the BGK kinetic model, which 
has a simple explicit solution in the steady state 0, |(| , 
does a very good job at predicting the transport proper- 
ties, especially in the case of the EHS system. An excep- 



tion is the normal stress difference (P z 



yy,s 



)/nT s , 



which vanishes in the BGK model but takes on (small) 
positive values in the simulations. The good agreement 
between simulation data for the elements of the pressure 
tensor and the BGK predictions was already noted in 
Ref. [2(|, although there the kinetic model was slightly 
different from the one considered here. 

Figure ^] also includes results obtained from control 
simulations carried out in the case a = 0.5 on EHS but 
with a refined friction constant 7 = |Co(^i2)/(^i2)o in- 
stead of the simple one 7 = ^Co- We observe that the 
agreement with the IHS data improves for those quanti- 
ties that were already reasonably well described by the 
simple EHS system, namely T s , P xy , s , and P xx ^ s - P yy , s - 
On the other hand, the delicate normal stress difference 



Pz 



which otherwise is quite small, is still about 



40% smaller in both EHS systems than in the IHS system. 
This indicates that whenever the discrepancies between 
the IHS and EHS results are relatively important, they 
are hardly affected by a more sophisticated choice of the 
friction constant 7. 

From a rheological point of view, it is worthwhile intro- 
ducing the nonlinear shear viscosity 77* = — (P xy /a)/r)o = 
— (P X y/nT)/a* and the viscometric functions \&i = 
(P yy - P xx )/nTa* 2 and f 2 = (P zz - P yy )/nTa* 2 . In 
the hydrodynamic transient regime, they are functions 
of a* (for a given value of a), as was illustrated in Fig. 
ITU1 in the case of 77* at a = 0.5. In the steady state 
(a* — > a*) these quantities become functions of a only. 
Equivalently, by eliminating a in favor of the steady-state 
reduced shear rate a* , those rheological quantities can be 
seen as functions of a* . This is the representation shown 
in Fig. El We observe that the curve rj*(a*) for EHS 
is very close to the one for IHS, even though the EHS 
points are shifted with respect to the IHS points corre- 
sponding to the same value of a, the shift increasing as 
a decreases. On the other hand, the viscometric effects 
are more pronounced in the IHS case than in the EHS 
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FIG. 15: (Color online) Steady-state values of the reduced 
shear viscosity 77* = 77/770 and the viscometric functions 
-*i = (Pxx - P yy )/nTa* 2 and * 2 = {P,z - P yy )/nTa* 2 
as functions of the reduced shear rate a^. The open symbols 
are simulation results for IHS and EHS and two values of the 
shear rate, while the dotted lines correspond to the solution 
of a BGK model. Note that in the latter model *&2,s = 0. 
The filled diamonds represent simulation data of EHS with 
a friction constant 7 = ^olV^) / ^12)0 in the case a = 0.5, 
ar° = 4. 



case. Paradoxically, the BGK curves for 77* and \&i, s are 
generally closer to IHS than to EHS in the representation 
of Fig. E| It is also noteworthy that the isolated points 
corresponding to EHS with 7 = |Co(^i2)/(^i2)o seems 
to be consistent with the curve obtained by joining the 
other EHS points. 

The cv-dependence of (Vi 2 ) / (Vi 2 )o and the cumulants 
C12 and 03 in the steady state is displayed in Fig. ^j] We 
recall that in the IHS system the ratio (Vi 2 ) / (Vi 2 )o co- 
incides with the ratio C/Co between the true cooling rate 
and its local equilibrium estimate. It is observed that the 
local equilibrium approximation underestimates the cool- 
ing rate by a few percent, essentially due to the distortion 
induced by the shearing. This distortion is well captured 
by the EHS system, despite the fact that its cooling rate 
is, by construction, given by Co- As said above, the posi- 
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FIG. 16: (Color online) Steady-state values of the the ratio 
(^i 3 2)/(^i2)o, the fourth cumulant 0,2, and the sixth cumulant 
—03 as functions of the coefficient of restitution a. The open 
symbols are simulation results for IHS and EHS and two val- 
ues of the shear rate, while the dotted lines in the middle and 
bottom panels correspond to the solution of a BGK model. 
The filled diamonds represent simulation data of EHS with 
a friction constant 7 = |Co(Vi 3 2)/{Vi 3 2 )o in the case a = 0.5, 



tive values of a 2 and, especially, —03 are indicators of an 
overpopulation effect (with respect to the Maxwellian) of 
the high- velocity tail of the distribution, this effect being 
stronger than and practically independent of the typical 
high- velocity overpopulation of granular gases in uniform 
and isotropic states |U|53- As a matter of fact, the dis- 
tribution function of the (frictional) EHS gas in the ho- 
mogeneous cooling state, as well as in the state heated by 
a white-noise forcing, is exactly a Gaussian. Therefore, 
the overpopulation measured by a 2 and —03 is basically 
a shearing effect. The fourth cumulant is well accounted 
for by the EHS gas, but the magnitude of the sixth cumu- 
lant is larger for IHS than for EHS. As happens with the 
normal stress difference P zz — P yy (see the bottom panel 
of Fig. 1 1 41) . the cumulant 0,3 is a very sensitive quantity 
that probes subtle details of the IHS velocity distribu- 
tion function not sufficiently well captured by the EHS 
system for a < 0.85. 



It is worthwhile noting that (V^)/{V^)o and a 2 are 
slightly larger in the EHS case than in the IHS case 
for a > 0.7, while the opposite happens for a < 0.7. 
This is reminiscent of the situation in the homogeneous 
cooling state and in the white- noise heated state |5lj . 
where a 2 HS < af HS = for a > V2/2 ~ 0.71, while 
4 HS > of HS = for a < V2/2. In fact, we have observed 
that the difference a 2 HS — af HS in the USF is quite close 
to the difference in the white-noise heated state. Figure 
Hoi also shows a strong correlation between the values of 
(V12)/ (^12)0 and those of 02- More precisely, we have 
empirically checked that (Vi 2 )/{Vi 2 )o — 1 + 0.27o 2 both 
for IHS and EHS in USF, in contrast to (V? 2 )/(Vf 2 ) ~ 
1 + jqO, 2 in isotropic states [5j . 

In Fig. E| we have included the BGK predictions for 
the cumulants a 2 and 03. Not being associated with con- 
ventional velocity moments, the evaluation of (V^) [cf. 
Eq. I|2.11(l ] from the BGK distribution function requires 
heavy numerical work and so it is not included in Fig. 1161 
We observe that the increase of a 2 with increasing inelas- 
ticity is well captured by the BGK model. On the other 
hand, the a-dependence of —03 is described by the BGK 
model at a qualitative level only, predicting in general 
too high values. 

Regarding the control simulations for EHS with a — 
0.5 and 7 = |Co{^i2)/(^i2)o; we observe that the values 
°f (^i2)/(^i2)o and a 2 are not actually improved, while 
the improvement on 0,3 is very small. 

The quantities plotted in Figs. I14H10I provide use- 
ful (indirect) information about the velocity distribution 
function of the steady-state USF. However, they essen- 
tially probe the domain of low and moderate velocities 
(say V < 2y/2T/m), except perhaps in the cases of 
a 2 and, especially, 03, which are more sensitive to the 
high-velocity tail. In order to analyze more directly the 
shape of the velocity distribution, its anisotropy, and 
the high-velocity tail, we have measured in the simula- 
tions the steady-state marginal distributions defined by 
Eqs. H3.18fl - I|3.20|) . As representative examples, Figs.lT7l 
(linear scale) and 1181 (logarithmic scale) show <j4 ((7 X ), 

9y + \Cy)i and F(C) for a = 0.9 and a = .5. We have 
checked that the symmetry properties l|3.21|l arc fulfilled 
and that the curves obtained from the two values of the 
shear rate (ar° = 4 and ar° = 0.1) practically coincide. 
In fact, to improve the statistics, the simulation data rep- 
resented in Fig. E] have been averaged over both shear 
rates and, in addition, the symmetry properties (|3.21|) 
have been exploited to make 
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Two anisotropic features of the USF state are quite ap- 



parent. First, the functions gi + \C x ) and g y ^'{C y ) are 
clearly asymmetric, namely ffa (|Ca.|) < g x + \~\C x \) 



,(+)/ 



and g y + \\C y \) < g y +> (-\C y \). This is a physical ef- 
fect induced by the shearing, in consistency with P xy oc 
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FIG. 17: (Color online) Linear plots of the marginal velocity distribution functions gi + \C x ), gy + \C y ), and F(C) for a — 0.9 
(left panels) and a — 0.5 (right panels). The solid and dashed lines represent simulation results for IHS and EHS, respectively, 
the dash-dot lines are the BGK predictions, and the dotted lines are the (local) equilibrium distributions. 
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distribution gi +) is broader than g y , in agreement with 
the fact that P xx - P m oc (C*) - (C?) > 0. These two 
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effects are obviously more pronounced for a 
for a = 0.9. 



Figure El shows that an excellent agreement between 
the IHS and EHS distributions for an inelasticity a = 0.9 
exits in the region of low and intermediate velocities, 
in consistency with the results displayed in Figs. 1141 - 



ITSl For a = 0.5, it is observed that the EHS distri- 
butions g x + \c x ), g^~\C y ), and F(C) are more popu- 
lated than the IHS ones in the regions — 1 < C x < 0, 
—0.3 < C y < 0, and C < 0.5, respectively. It is also 
interesting to note that, in the region of thermal veloci- 
ties, the orientation-averaged distribution function F(C) 
is much less distorted with respect to the Maxwellian 
than the marginal distributions g x + \c x ) and gy + \c y ), 
especially in the case a = 0.9. Another important find- 
ing is that the BGK kinetic model, not only captures well 
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FIG. 18: (Color online) Logarithmic plots of the marginal velocity distribution functions ffi (C x ), <?y (C H ), and F(C) for 
a = 0.9 (left panels) and a = 0.5 (right panels). The solid and dashed lines represent simulation results for IHS and EHS, 
respectively, the dash-dot lines are the BGK predictions, and the dotted lines are the (local) equilibrium distributions. 



"global" or average properties, such as the hydrodynamic 
quantities (cf. Figs. I14H16|I . but also the "local" details 
of the velocity distribution function. 

The logarithmic scale employed in Fig. ED has been 
chosen to reveal the high-velocity tails of the distribu- 
tions. We observe an overpopulation of both tails of 
gi + \ especially at a = 0.5. In the case of gy + \ how- 
ever, the overpopulation seems to affect the tail C y < 
only. The high-energy overpopulation is clearly apparent 
in the distribution of the magnitude of the velocity, F(C). 
In a recent paper, Bobylev et al. |53| have proven that 



In F(C) ~ — C M , with fi > 1, for asymptotically large ve- 
locities in the USF. Although it is not clear whether our 
simulation data have reached the high-velocity regime 
where the asymptotic law lnF(C) ~ — C p dominates, 
the bottom panels of Fig. ED seem to be consistent with 
this law with fi ~ 1. 

The comments in the preceding paragraph apply 
equally to IHS and EHS. As a matter of fact, the sim- 
ulation data for both systems in the case a = 0.9 are 
hardly distinguishable for that value of a. In the case 
a = 0.5, however, the high-velocity values of <)4 + (C x ) 



20 



and F(C) are larger for IHS than for EHS, while the op- 
posite happens for the high- velocity values of gy + \C y ). 
As for the BGK distribution function, it cannot be ex- 
pected to be accurate beyond the domain of thermal ve- 
locities. This is confirmed by Fig. ^1 where we can ob- 
serve that at a = 0.9 the BGK model strongly exagger- 
ates the high- velocity overpopulation effects of gx(C x ) 
and F(C). This produces too large a value of the sixth 
cumulant —03, as observed in Fig. 1161 On the other hand, 
at a = 0.5 the BGK value of — a 3 agrees by accident 
with that of the IHS system and this explains why in 
that case the overpopulated tail predicted by the BGK 
model agrees well with that of the IHS system. Never- 
theless, the tail of g£ (C x ) for C x > as well as that of 
g y + \C y ) for C y < are strongly underestimated by the 
BGK model at a = 0.5. 

For the sake of clarity, we have not included in the 
right panels of Figs. 1171 and!18l the curves corresponding 
to our simulations of the EHS system with a = 0.5 and 
a friction constant 7 = ^Co(U 1 3 2)/(U 1 3 2 )o. In any case, 
the results are very close to those corresponding to the 
friction constant 7 = i(o- This shows that the main 
quantitative differences between highly dissipative IHS 
and EHS systems, in particular the high- velocity tails, 
are intrinsic to the different dynamics of both systems 
and so they are not avoided by a fine-tuning of the drag 
force acting on the elastic particles. 



VI. CONCLUSIONS 

In this paper we have dealt with the main nonequilib- 
rium properties of two classes of dissipative gases sub- 
ject to the so-called simple or uniform shear flow (USF). 
In the first class, the system is made of inelastic hard 
spheres (IHS) with a constant coefficient of restitution 
a < 1 . The inelasticity of collisions provides an internal 
energy sink, characterized by a cooling rate £. In the 
second class, the particles are elastic hard spheres (EHS) 
under the action of a drag force, so that the energy sink is 
now external and characterized by a friction coefficient 7. 
At a given value of a, some of the basic properties of the 
Boltzmann equation for the IHS gas are formally similar 
to those for the (frictional) EHS gas (4] if (i) the friction 
coefficient is chosen as a function of the local density and 
temperature, 7 = |Co x nT x l 2 {\ — a 2 ), and (ii) the col- 
lision rate of the EHS gas is /3(a) = s(l + a) times the 
collision rate of the IHS gas under the same conditions. 
The boundary conditions of the USF state provides an 
energy source (viscous heating) that competes with the 
energy dissipation (either collisional or frictional), until 
a steady state is eventually reached, resulting from the 
balance between both effects. 

The two main points we have intended to address in 
this paper are the following ones. On the one hand, we 
wanted to perform an extensive study of the physical 
properties of dissipative gases under USF, for both the 



transient and the steady states. While the USF state has 
been widely considered in the literature of granular flu- 
ids, some of its relevant aspects (hydrodynamic transient 
stage, cumulants, high-velocity tail, . . . ) have received 
little attention. As a second point, taking the USF as 
a paradigmatic and conceptually simple nonequilibrium 
state, we were interested in elucidating to what extent the 
EHS gas mimics the physical properties of the genuine 
IHS gas. To meet these goals, we have carried out com- 
puter simulations by the DSMC method on both classes 
of systems for ten values of the coefficient of restitution 
(a = 0.5-0.95 with a step Aa — 0.05) and two values 
of the shear rate (a — 4/t° and a = 0.1/r°). Below we 
summarize the main conclusions derived from our study. 

The duration of the transient period, when measured 
by the number of collisions per particle (s), is hardly de- 
pendent on the imposed shear rate a or on the initial 
state, but strongly depends on the coefficient of restitu- 
tion a. The larger the dissipation, the smaller the num- 
ber of collisions needed to reach the steady state. For 
instance, s ~ 40 at a = 0.9, while s ~ 20 at a = 0.5. 
Nevertheless, when the duration is measured by an exter- 
nal clock (for instance, in units of the initial mean free 
time), it becomes weakly dependent on a but strongly 
dependent on a: the smaller the imposed shear rate, the 
longer the transient period. 

The evolution toward the steady state proceeds in two 
stages. The first (kinetic) stage depends heavily on the 
initial preparation of the system and lasts a few colli- 
sions per particle. This is followed by a much slower 
hydrodynamic regime that becomes less and less depen- 
dent on the initial state. Once conveniently scaled with 
the thermal speed Vo(t) = y/2T(t)/m, the velocity dis- 
tribution function in the hydrodynamic regime depends 
on time through the reduced velocity C(t) = ~V/vo(t) 
and the reduced shear rate a*(t) cx a /y/TJfy only 0. In 
particular, at a given value of a, the (reduced) nonlin- 
ear shear viscosity rj* (a*) moves on a certain rheological 
curve, the steady-state value rj* = rj*(a*) representing 
just one point. This point splits the curve rj*(a*) into 
two branches. The branch a* < a* is accessible from 
initial states such that the dissipative cooling dominates 
over the viscous heating, so that T(t) decreases and a*(t) 
increases during the transient period. Conversely, the 
branch a* > a* corresponds to initial states where the 
viscous heating dominates, so that T(t) increases and 
a*(t) decreases, until the steady state is reached. 

The actual number of collisions per particle s(t) is 
slightly smaller than the local equilibrium estimate so(t). 
This indicates that (V) < (V) , whereas (V 2 ) = (V 2 ) = 
3T(t)/m by definition. The inequality (V) < (V)q is con- 
sistent with an underpopulation (with respect to the lo- 
cal equilibrium distribution) for moderate velocities that 
must be compensated by a high-velocity overpopulation. 
This qualitative reasoning is confirmed by the inequali- 
ties (V&) > (ViDo, (V 4 ) > (U 4 ) , and (V 6 ) > (U 6 ) ob- 
served in the simulations. In addition to these effects, the 
shearing motion induces a strong anisotropy in the nor- 
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mal stresses, namely P xx {t) > nT{t) > P zz (t) > P yy (t). 
In other words, a breakdown of the energy equipartition 
occurs, whereby the "temperature" associated with the 
degree of freedom parallel to the fluid motion is signif- 
icantly larger than that associated with the other two 
degrees of freedom as alr eady observed inprevious stud- 
ies pi e e m n m 0, m m miH m. 

We have paid special attention to the properties of the 
steady state, which is intrinsically independent of the im- 
posed shear rate and of the initial state. As expected, 
the distortion from the local equilibrium state (as mea- 
sured by the shear stress, the normal stress differences, 
the cumulants, . . . ) increases with the dissipation. This 
distortion is made quite apparent by the shapes of the 
steady-state (marginal) velocity distributions defined by 
Eqs. (|3.18l) - (|3.20|) . In particular, the high- velocity tails 

of s4 (C x ) and F{C) seem to be consistent with an ex- 
ponential overpopulation. 

All the above comments apply equally well to the IHS 
and EHS gases. Therefore, the main nonequilibrium and 
transport properties of a true IHS gas in the USF state 
are satisfactorily mimicked by an "equivalent" EHS gas. 
If one focus on the basic properties (say the steady-state 
reduced shear rate a*, shear stress P xy>a /nT s , or second 
cumulant a 2 ) it is almost impossible to distinguish the 
EHS values from the IHS ones if a > 0.7, the differences 
still being relatively small if a < 0.7. We have observed 
that {V12) / {Vi 2 )o and a 2 are in the EHS system slightly 
larger than in the IHS system if a > 0.7, while the oppo- 
site happens if a < 0.7; this is analogous to what happens 
in the homogeneous cooling state and in the white-noise 
heated state, in which cases the EHS distribution is ex- 
actly a Maxwellian. It is interesting to note that, even 
at a — 0.5, the full hydrodynamic curve rj*(a*) coin- 
cides for both systems, except that the location of the 
steady-state point 77* = rj*(a*) slightly changes. More 
delicate quantities (such as the normal stress differences 
or the sixth cumulant) keep being practically the same in 
both systems if a > 0.85. Even at the level of the veloc- 
ity distribution function itself, the EHS and IHS curves 
practically overlap (at least for the domain of velocities 
C < 6 accessible to our computer simulations) at a coef- 
ficient of restitution as realistic as a — 0.9. At a = 0.5, 
however, the distribution functions g x + \C' x ) and F(C) 
of the IHS system exhibit a visibly larger high-velocity 
overpopulation than those of the EHS system. 

As said above, in the EHS systems we have chosen the 
friction coefficient as 7 = i(o 0(1 «T 1 / 2 (1 — a 2 ), so that it 
is a functional of the distribution function only through 
the local density and temperature and, moreover, its de- 
pendence on a is explicit. Given that the true cooling 
rate is slightly larger than the local equilibrium estimate, 
i.e., C/Co = (^i2)/(^i2)o > 1) the imitation of the inelas- 
tic cooling rate by the EHS is not perfect. Therefore, one 
might reasonably expect that the discrepancies between 
the EHS and IHS results would diminish if the friction co- 
efficient were taken as 7 = §£0(^12)/ (^12)0- To test this 
expectation, we have performed complementary simula- 



tions of the EHS system with this more refined value of 
7 in the case of highest dissipation, i.e., a = 0.5. The re- 
sults show that, whenever the former agreement between 
EHS and IHS was fair, the new agreement is generally 
even better. However, those quantities (such as the nor- 
mal stress difference P zz — P yy and the sixth cumulant 
03) that turned out to be especially sensitive to the dissi- 
pation mechanism (collisional inelasticity versus external 
friction) are practically unaffected by the new choice of 
7. Moreover, the high-velocity tails of both versions of 
the EHS gas are practically the same, being both smaller 
than the IHS tail. Therefore, the (subtle) discrepancies 
between the IHS and EHS systems in cases of high dissi- 
pation (say a < 0.7) seem to be intrinsic to their distinct 
dynamics. Taking this into account, there is no practical 
reason to propose for the drag force acting on the EHS 
a friction coefficient different from the local equilibrium 
value 7 = i£o- Concerning the collision rate coefficient 
(3(a), the choice (3{a) = |(1 +a) is recommended by cri- 
teria of simplicity and consistency with the cases of mix- 
tures and dense gases Q . Moreover, we have checked (not 
shown in this paper) that an alternative choice, namely 
(3{a) — i(l + a)(2 + a), although reproducing well the 
Navier-Stokes shear viscosity Q , provides results for the 
nonlinear shear viscosity in worse agreement with the IHS 
ones than those reported here with f3(a) — |(1 + a). 

The (approximate) equivalence IHS^->EHS can be used 
to transfer to granular gases part of the expertise accumu- 
lated for a long time on the kinetic theory of elastic par- 
ticles. In particular, the celebrated Bhatnagar-Gross- 
Krook (BGK) kinetic model of the Boltzmann equation 
can be readily extended to granular gases 0, • In this 
paper we have compared the solution of the BGK model 
for USF H H with the simulation data. While it is 
generally believed that the BGK model would be accu- 
rate only for states near equilibrium and/or in the quasi- 
elastic limit, our results show that, despite its simplic- 
ity, the model succeeds in capturing quantitatively the 
evolution and steady-state values of the main transport 
properties (temperature, shear stress, and normal stress 
difference P xx — P yy ) and even of the fourth cumulant a 2 . 
However, the small nonzero difference P zz — P yy is not 
accounted for by the BGK model and the sixth cumu- 
lant 03 agrees with the simulation data at a qualitative 
level only. All of this is consistent with the observation 
that the BGK velocity distribution function is reliable in 
the thermal region (C < 2), but not in the high- velocity 
domain. 

In summary, we conclude that the solutions of the 
Boltzmann equations for IHS are very similar to those 
for EHS (in the latter case with a smaller collision rate 
and under the action of an adequate drag force). Thus 
the temporal evolution toward the steady state and the 
properties of the latter are mainly governed by the com- 
mon feature of energy dissipation, without any significant 
influence of the detailed mechanism behind it. Only for 
very high dissipation (say a < 0.7) and for properties 
probing the velocity domain beyond the thermal region, 
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does the IHS system imprint its signature and distinguish 
from the "disguised" EHS system. In this paper we have 
restricted ourselves to the USF, but we plan to perform 
a similar comparison in other states, especially in those 
where the heat flux, rather than the pressure tensor, is 
the relevant quantity. We will also undertake a parallel 
study in the case of dilute mixtures, as well as in the 
case of dense gases (complemented by molecular dynam- 
ics simulations), following the schemes discussed in the 
preceding paper Q. 
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